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Summary 


This final report will document the accomplishments of the work of this project. 

1. The incremental-iterative (II) form of the reverse-mode (adjoint) method for computing first- 
order (FO) aerodynamic sensitivity derivatives (SDs) has been successfully implemented 
and tested in a 2D CFD code (called ANSERS) using the reverse-mode capability of 
ADIFOR 3.0. These preceding results compared very well with similar SDs computed via a 
black-box (BB) application of the reverse-mode capability of ADIFOR 3.0, and also with 
similar SDs calculated via the method of finite differences. Detailed documentation of all of 
the preceding is provided in Appendix A. 

2. Second-order (SO) SDs have been implemented in the 2D ASNWERS code using the very 
efficient strategy that was originally proposed (but not previously tested) of Reference 3, 
Appendix A. Furthermore, these SO SOs have been validated for accuracy and computational 
efficiency; detailed documentation of these SO SDs is provided in Appendix A. 

3. Studies were conducted in Quasi- ID and 2D concerning the “smoothness” (or lack of 
smoothness) of the FO and SO SD’s for flows with shock waves. The phenomenon is 
documented in the publications of this study (listed subsequently), however, the specific 
numerical mechanism which is responsible for this unsmoothness phenomenon was not 
discovered. 

4. The FO and SO derivatives for Quasi- ID and 2D flows were applied to predict aerodynamic 
design uncertainties, and were also applied in robust design optimization studies. Detailed 
Documentation of this is provided in Appendix B. 
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An *ger«egta M ter a ttv g gnnmacti differentiating advanced flow codes is successfully demonstrated 

on a twtHt iB W Btom i tevfedd nodd problem. The method employs tbe reverse-mode capabflfty of tbe tatooutk- 
differentiation software tool ADIFOR 3.0 and is proves to yield accurate first-order aerodynamic sensitivity 
derivatives. A substantial reduction in CPU time aad computer memory Is demonstrated fas comparison with 
results from a stra i g htfor ward, black-box reverse-mode application of ADIFOR 3.0 to tbe same flow code. An 
ADIFQR-assfatcd procedure for accurate second-order aerodynamic sensitivity derivatives is successful verified 
on an ioviscid transonic lilting airfoil example problem. The method requires that first-order derivatives are 
calculated first using both the forward (direct) and reverse (adjoint) procedures; then, a very efficient noniterative 
calculation of aO second-order derivatives can be accomplished. Accurate second derivatives (Le^ the complete 
Hessian matrices) of lift, wave drag, and pitching-moment coefficients are calculated with respect to geometric 
shape, angle of attack, and freestream Mach number. 


L Introduction 

T HIS paper revisits and focuses entirely on the computational 
challenges that are associated with the efficient calculation of 
aerodynamic sensitivity derivatives (SDs) from advanced computa- 
tional fluid dynamics (CFD) codes. Of course, an accurate efficient 
methodology for obtaining these derivatives is a critical prerequi- 
site concern that must be addressed first by the aerodynamic design 
engineer who chooses any gradient- based method(s) for design op- 
timization and/or for estimating quantities related to aerodynamic 
uncertainty. Thus computing SDs from high-fidelity nonlinear CFD 
codes is an enabling technology for design of advanced concept 
vehicles. 

In recent years significant progress has been achieved in the ef- 
ficient calculation of accurate SDs from these CFD codes. 1 The 
automatic differentiation (AD) software tool ADIFOR (Automatic 
Differentiation of FORTRAN) has been proven an effective tod for 
extracting aerodynamic SDs from these modem CFD codes. 2-6 The 
foundation of the present work is found in Refs. 3 and 6; the present 
study builds on these earlier studies in an effort to exploit the full 
potential of the latest version of ADIFOR 3.0 (Ref. 7) for obtaining 
SDs from CFD codes. 

In Ref. 2 a strategy known as the ADD method was first pro- 
posed and later successfully demonstrated in Ref. 3, whereby AD 
was applied to a CFD code in incremental-iterative (I-I) form. The 
ADD method is a hybrid (compromise) scheme, designed to main- 
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tain as much as possible the computational efficiency of a hand- 
differentiated (HD) approach and the ease of implementation of 
a straightforward black-box (BB) application of AD; at the same 
time the accuracy of the SDs is not compromised. A comprehensive 
overview of the development of the ADD scheme is given in Ref. 
3. Also included in Ref. 3 is a comparison of the ADD scheme with 
the HD and BB approaches; computational issues associated with 
CPU time, computer memory, and SD accuracy are discussed. The 
two-dimensional effort of Ref. 3 was later extended to the three- 
dimensional code CFL3D, including “in-parallel” computation of 
the derivatives. 8 * 9 Appropriate references to the version of CFL3D 
used can be found in Ref. 6. 

The success reported in these previous works 3 ’ 8,9 could be con- 
sidered limited, however, because all ADIFOR implementations 
reported therein were forward-mode (direct) differentiations. It is 
very difficult to make any forward-mode implementation of deriva- 
tive calculations computationally competitive with a reverse-mode 
(adjoint) implementation whenever the number of design variables 
(NDV) of interest is considerably larger than the number of out- 
put functions (NOF) of interest, and NDV much greater than NOF 
is more typical for aerodynamic design problems. In recent stud- 
ies the new reverse-mode capability of ADIFOR 3.0 (not available 
for tire earlier referenced studies) has been successfufiy verified in 
Ref. 6 by application to a parallel version of CFL3D and in Ref. 
10 by application to a sequential linear aerodynamics code. These 
applications resulted in accurate design SDs as well as stability and 
control derivatives, respectively. The application reported in Ref. 6 
involved BB AD of the entire CFD code, but iterative execution of 
the reverse mode was required only over the last iteration of the 
function evaluation. 

In the present study it is proposed and demonstrated that the 
reverse-mode capability of ADIFOR 3.0 can also be applied to CFD 
codes in I-I form, resulting in a hybrid adjoint-variable ( AV) scheme 
(known herein as the ADII-AV method) that is analogous to the 
forward-mode ADD scheme of Ref. 3 and elsewhere. The motiva- 
tion of this new reverse-mode ADII-AV scheme is identical to that 
of the earlier forward-mode ADD method: greater computational 
efficiency is sought over a BB implementation of AD, without any 
loss of accuracy in the calculated SDs and without unmanageable 
complications upon implementation. 
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Following development of the proposed new ADH-AV scheme, 
the second focus of the present study is that of calculating second- 
order (SO) aerodynamic SDs from CFD codes. The motivation for 
calculating SO SDs is to advance the possibility of (or greater capa- 
bility of) 1) second-order gradient-based aerodynamic design op- 
timization, 2) analysis and design involving vehicle stability and 
control and 3) robust design [i.e., design under uncertainty, where 
a first-order, second-moment method requires SO SDs 11 ]. This sec- 
ond part of the present study is another extension of Ref. 3, wherein 
the computational issues associated with calculating these higher- 
order derivatives were addressed, and sample calculations of SO 
derivatives using AD were repented from a two-dimensional CFD 
code. 

In Ref. 3, four procedures few calculating SO CFD SDs were pro- 
posed, but only one of the less efficient methods was actually tested; 
ADIFOR 3.0 currently provides three forward-mode variations for 
the calculation of SO SDs by similarly inefficient methods. The most 
efficient (for large NDV) SO SD scheme was not tested in the earlier 
study 3 but has been successfully implemented in the present study. 
Reverse-mode (adjoint-based) differentiation is required within this 
efficient SO SD scheme, via either HD or AD. However, with the 
availability of ADIFOR 3.0 and the new ADII-AV scheme the door 
has been opened few AD implementation and testing of this SO SD 
scheme few CFD codes. The results of this effort to date are reported 
here. These efficiently computed SO SDs have been used to demon- 
strate an approach for CFD input uncertainty propagation and robust 
design optimization for a quasi-one-dimensional flow application in 
Refill. 

IL Basic Equations and Theoretical Development 

The equations summarized subsequently are discussed in greater 
detail in the references, in particular, Refi 3. These concepts are 
known in the mathematical optimization community, 12 but the de- 
tails developed here do not appear to be generally known throughout 
the CFD community. The aerodynamic output functions of interest 
F and the discretized conservation laws of steady compressible fluid 
flow If, including boundary conditions, can be represented symbol- 
ically as follows. 

Aerodynamic output functions: 


In the preceding equations i, j 9 and / are “free” indices, and re- 
peated indices m and p are (by convention) “summation” indices. 
The reverse-mode (adjoint) approach for the FO SDs is developed 
starting with an application of the chain rule: 


dFi dR t _ dFt 
dRi dQ m ~ dQ m 


in order to show that the conventional adjoint variable is indeed a 
derivative. With Eq. (5) it then follows from Eq. (4) that Eq. (3) can 
be written as 



M±(™L Y . j-El 

dR, (a*, pj dbj) dX„ pi db 


j 


A more conventional derivation of the adjoint-variable method (i.e., 
the Lagrange-multiplier method) gives 


dFi _ . / 3R, , SR'\ dFi dFi 

F " = ~ W ^6“ (7) 


dR, dFi _ fl 

Gim= ll W~ + W n ~ in 


( 8 ) 


where A.,/ is called the adjoint variable. Comparison of Eqs. (6) and 
(5) with Eqs. (7) and (8), respectively, reveals the identity 


dRi 




(9) 


One objective of this particular development of the AV method for 
aerodynamic SDs is to ensure that the relationship given by Eq. (9) 
is dearly understood; that is, X u is the derivative of the output F, 
with respect to the intermediate variable /?/ and is accumulated in 
the reverse-mode AD. 

The I-I strategies for solving the preceding equations for Q f mJ 
and/or X a required for F[j are reviewed here; additional detail is 
found in Refs. 3 and 13 and elsewhere. The I-I method for solving 
the nonlinear flow of Eq. (2) is 


F = F[g(*), *(*),&] (1) 


e" +1 = ei:-pjX 


( 10 ) 


Nonlinear state equations: 

/f=g[e(«,m« = o (2) 

where Q is the vector of state (field) variables, X is the vector of 
computational grid coordinates, and Ms the vector of input (design) 
variables. 

A. First-Order Sensitivity Derivatives 

Subject to the following definitions, index (summation) notation 
is now introduced: 



(This notation will be necessary to avoid subsequent ambiguity when 
the SO SD methods are presented.) The forward-mode (direct) ap- 
proach for calculating first-order (FO) SDs is developed by differ- 
entiation of Eqs. ( 1 ) and (2) with respect to die design variables; the 
result is 


F' — 


dFj 

dbj 



3F, dFi 

Wj, pi + Wj 



dR, 

3Qm 


Qmj + 


dR, 

dX p 


yf . 8R l _ 0 

Xp > + W i ~ 11 


(3) 

(4) 


where the superscript N is the iteration (pseudo-time step) index 
and the operator 


p n ae* 

^ dR f 



(ID 


represents tire solution algorithm of the particular CFD code of 
choice . The tilde in Eq. (11) serves to indicate that can be viewed 

as any computationally efficient approximation (often a very crude 
approximation) of the exact operator associated with true Newton- 
Raphson iteration. Thus the CFD solution algorithm is simply quasi- 
Newton iteration. 

The I-I method for solving the forward mode FO SD equation (4) 
is 


Q 


fH + l 

mj 



-PmlR 


U 


where superscript M is the FO SD iteration index and 


R 


,M _ dRl 
" “ 3Q m 


Qmj 


+ 


3R,_ Y , dR, 
dX„ pJ + dbj 


( 12 ) 


(13) 


With the I-I methodology the CFD flow solution operator F„i is 
also used to solve the SD equations; this operator in Eq. (12) is 
evaluated and fixed using the steady-state solution for the nonlinear 
flow. The requisite terms of Eq. (13) are constructed either by hand 
differentiation (i.e., the HDII method, which is very tedious and time 
consuming to complete with accuracy for advanced CFD codes) or 
by AD, which is the forward-mode ADII method of previous studies. 
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In contrast with the ADQ method, a straightforward BB appli- 
cation of AD to the CFD code, which is the ADBB method, is 
represented symbolically as 


' = ^ = W r , _ _ /dG\ 

d»* \&) llk 3M*/' “ d*t 


C 1 = C «/ - O v - KjX? (14) 

Oeariy ADII [Eqs. (12) and (13)] and ADBB [Eq. (14)] yield the 
same result at steady-state convergence of each [recall Eq. (2)]; 
however, ADII is potentially more efficient than ADBB because of 
user intervention in the application of AD. With ADII the following 
is true: 

1) The operator can be evaluated only once [hence denoted 
Pmi in Eq. (12)] using the steady-state field variables Q and then 
reused for all M iterations and for all j = NDV design variables in 
obtaining the Q f mJ , 

2) All derivatives except Q' mj can be computed once outride the 
iteration loop and frozen far reuse inride the loop. 

3) Evaluation of die terms P'^j in Eq. (14) can be avoided com- 
pletely for all iterations and all design variables. 

The I-I method for solving the reverse-mode AV, FO SD Eq. (9) 
is 


where 


p r M 

*il — k il ~ 


* BRj BFj 
" 32m 32m 


(15) 

(16) 


The requisite terms of Eq. (16) arc constructed either by hand (Le., 
the HDII-AV method, having the same drawbacks as the forward- 
mode HDII method) or by AD, which is the proposed new ADD- AV 
scheme. The BB AD in reverse mode (the ADBB-AV method) has 
been verified in Ref. 6. The objective of the proposed ADD-AV 
scheme is improved computational efficiency over the ADBB-AV 
approach without resulting loss of accuracy or significant loss in 
the ease of implementation. The mechanisms from which improved 
computational efficiency can be expected are analogous to those 
explained before when the forward-mode ADD and ADBB methods 
were contrasted. Furthermore, the ADII-AV scheme should lend 
itself to more permanent generalized coding implementations than 
the ADBB-AV approach. This is because with the ADD-AV method 
the manner in which AD is applied is independent of, yet valid for, 
all of the particular aerodynamic inputs and outputs of interest 
The forward-mode application of ADIFOR produces FORTRAN 
source code for very efficient calculation of the vector (or ma- 
trix) product that results from the postmultiplication of a large 
Jacobian matrix by a known input vector (or matrix). This at- 
tribute of forward-mode AD is exactly what was required to con- 
struct the ADII method; specifically, the terms {BRi/SQ m )Q f mj an d 
(dRi/BX p )X f P j of Eq. (13) are of this type. In contrast, however, 
the forward-mode application of ADIFOR produces source code 
that is prohibitively inefficient for calculating the premultiplication 
of a large Jacobian matrix by a known input vector (or matrix). This 
weakness of the forward-mode application of ADIFOR is exactly 
the strength of the reverse-mode option now available in ADIFOR 
3.0. Thus, the proposed new efficient ADII-AV scheme has become 
possible with this reverse-mode capability. That is, through reverse- 
mock application of ADIFOR 3.0 it is now feasible to produce (au- 
tomatically) the source code required for efficient evaluation of the 
term Xff ( BRi/BQ m ) in Eq. (16), that is, the prcmuitiplication of a 
large Jacobian matrix by a vector. 


B. Second-Order Sensitivity Derivatives 

The SO SD methods are presented in index notation subject to 


the following definitions: 


_ tfF, /<fF\ 

,Jk d b k dbj {<&)„' 

_ cfR, /d^X 

d htoj \ d^ 2 / ,j t 

„ _ c P2„ _(#Q\ 
mJk dfcd bj W) mJk 

_ d% (<?X\ 
ttk M k dbj \<&) fjk 


The following differential operator is also introduced for subsequent 
notational compactness: 


do = _30„, +®<>y» , £0 
D b k 32. dX 9 qk db k 


(17) 


where repeated indices n and q are summation indices. 

Differentiation of the FO forward-mode equations (3) and (4) 
with respect to the design variables yields SO method 1: 




tli - ELn" 

db k dbj 32™ Um/t ' lpjt 


+ ——X" 

i+ ax p "“ 


Db k 


(18) 


_ d 2 Ri 


dfy q 1 I dfy Y» 

32« ** »*/"* 


dr;, 

+ D6?=°"‘ (19 > 


The terms of DF'j/Dbt and D/^ /Db* are many and very com- 
plicated; detailed expansion of these terms is provided in the Ap- 
pendix. Using symmetry of the Hessian = Q'^j, SO method 1 

requires (NDV 2 + N DV)/2 solutions of the large linear systems 
of Eq. (19) for in addition, the method requires NDV solu- 
tions of Eq. (4) for the FO SDs Q ' mj . SO method 1 was verified for 
a two-dimensional CFD code in Ref. 3 by ADBB differentiation of 
the code’s existing HDII scheme [Eqs. (12) and (13)] for the FO 
SDs. 

Alternatively, differentiation of the FO reverse-mode, Eqs. (7-9), 
with respect to the design variables, yields SO method 2: 


d*F, /3R, 3R,\ 

tJk db k dbj llk \ax p Xl,J+ dbj) 


. (**i , , ^,\ D F!j 

+ {jr p +x “ar p ) x " k+ m 

(20) 

p/ dG,„ dRi DG im 

imk ~ d b k ~ "*32™ D b k ~ imt 

(21) 


SO method 2 requires NDV x NOF solutions of the large linear sys- 
tems of Eq. (21) for X\ lk \ in addition, the method requires NDV 
solutions of the FO equation (4) for Q' mj plus NOF solutions of the 
FO equation (9) for k it . This SO method 2 is eliminated from further 
consideration because it is unconditionally less computationally ef- 
ficient than the remaining two SO SD methods. 

Introduction of the AV approach within SO method 1 to eliminate 
Q'njk yields SO method 3: 


r . * a My* + 
Fi *~\aT p + Xi, aT p ) x ^ k + m 


+ Xu 


Db k 


( 22 ) 


SO method 4 is similar and computationally equivalent to SO 
method 3 and is developed by introduction of the AV approach 
within SO method 2 to eliminate X' ilk ; the result is the identity 


, _ n , DC*. 

" D b k Qmj D b k 


(23) 


where SO method 4 uses Eq. (23) to replace equivalent terms within 
SO method 3. The equivalent SO SD methods 3 and 4 do not re- 
quire solution of large systems of linear equations for higher-order 
derivatives such as Q n mjk or X' ilk . These two SO SD schemes do, 
however, require solution of both forward-mode and reverse-mode 
equations (4) and (8) for Q' mJ and X//, respectively. This is a total of 
only NDV +NOF solutions of large systems of linear equations. 

One significant conclusion of the preceding analysis is that SO 
methods 3 or 4 should be computationally more efficient when- 
ever NDV 2 +NDV is greater than 2 x NOF. With typical design 
problems in aerodynamics, NDV is often much larger than NOF; 
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typically, NOT is three or less (often only one), whereas NDV is on 
the order of tens to hundreds. The advantage in favor of method 3 
or 4 for SO SDs is then overwhelming because of the NDV 2 term, 
which dominates. Once both the forward-mode and reverse-mode 
schemes are in place for calculating the FO SDs, then complete 
SO SD information is available almost “for free,” that is, the SO 
SD are obtained through an explicit, noniterative calculation. The 
source code for implementation of method 3 or 4 is constructed 
“automatically” via BB application of the forward-mode capability 
of ADIFOR to appropriate pieces of the existing source code from 
which the FO SDs are obtained. For example, the extremely com- 
plex toms DFjj/DbtiVRy/Dbt, and/or DG im /Db k (see Appendix) 
of methods 3 and/or 4 are easily constructed with a forward-mode 
application of AD. 

In Ref. 3, SO methods 3 and 4 were proposed but not actually 
tested. Consequently, one primary goal erf the present study is suc- 
cessful implementation ami verification of the highly efficient SO 
method 3 (or equivalently, method 4); method 3 is actually chosen 
in this study. 

m. Results and Discussion 

A. First-Order Sensitivity Derivatives, 

ADII-AV Method, Model Problem 

The proposed ADII-AV method [Eqs. (15) and (16)] has been 
successfully implemented in a CFD code and verified for accuracy 
on a simple two-dimensional inviscid internal flow model problem. 
This CFD code solves the two-dimensional Euler equations by a 
conventional upwind finite volume approach on a very coarse grid 
but one that is sufficient for computationally verifying SDs. As ex- 
pected, when FO SDs computed by the new ADII-AV scheme are 
compared with SDs computed by a hand -differentiated implemen- 
tation of Eqs. (15) and (16) (i.e., the HDII-AV approach), the results 
are the same at convergence, as well as at each I-I step. In addition, 
the accuracy of the computed SDs has been successfully verified by 
a finite difference method. 

Preliminary timings were conducted on a Sun workstation to eval- 
uate the potential for improved computational efficiency of the new 
ADII-AV scheme with respect to the ADBB-AV approach of Ref. 6. 
Computational timing comparisons are given in Table 1 , which fo- 
cuses exclusively on AD performance. Therefore, relative timings 
are given as CPU time per iteration per grid point per differentiated- 
aerodynamic -output function. Furthermore, each timing result has 
been scaled by the comparable timing result obtained from the very 
efficient hand-differentiated reverse-mode scheme (i.e., the HDII- 
AV method). Table 1 illustrates that although the new ADII-AV 
scheme is almost five times slower than the efficient HDII-AV 
scheme it represents a substantial improvement over results obtained 
from the straightforward black-box procedure (i.e., ADBB-AV is 
about eight times slower than HDII-AV). 

The improvement in computational efficiency achieved to date 
is substantial when the reverse-mode application of ADIFOR 3.0 
in incremental-iterative form is compared with the black-box ap- 
proach. Furthermore, the timing result for the ADII-AV scheme is 
projected to improve by an additional 30% over that reported here. 
Thus the relative timing given in Tfcble I for ADD- AV/HDH-AV is 
projected to drop from 4.7 to about 3.3. This projection is based on 
using a strategy where the forward-pass execution of the ADIFOR- 
enhanced, reverse-mode code will be performed only once (instead 
of during each iteration) in order to create the required ADIFOR log 


Tablet Relative CPU timing 
comparison, model problem 


Reverse-mode 

Relative 

method tested 

timing* 

HDII-AV 

1.0 

ADII-AV 

4.7 

ADBB-AV 

7.9 


*CPU tkne/ilcMtio^g 
function. 


files. Thereafter, by repeatedly reusing these fixed log files oily re- 
verse passes will be repeatedly executed during the iterative solution 
process for all aerodynamic output functions of interest. 

Another important computational concern mitigated by the new 
ADII-AV method is computer memory, particularly the issue of large 
disk files created during execution of reverse-mode derivative code 
created by ADIFOR 3.0. With the black-box (ADBB-AV) approach 
these large ADIFOR log files (which are created on a forward-pass 
execution and are read during the reverse pass) will accumulate and 
become larger with every iteration of the ADIFOR-eohanced flow 
code. This file growth can rapidly deplete the available disk space, 
even on the largest computers. In Ref. 6 this difficulty was addressed 
by development of the iterated reverse-mode scheme, where only 
the log files for the final forward-pass iteration are stored and used 
during the subsequent iterative solution for the derivatives. With 
the ADII-AV approach, however, the required disk space is not as 
restrictive an issue because it remains fixed and does not accumulate 
during the iterative solution process. In the present example the total 
storage requirement for log files with the ADII-AV method is only 
64% of that required for a single iteration of the ADBB-AV method. 

In addition to the log-file disk memory, required only for the new 
reverse-mode capability of ADIFOR, there are substantial additional 
core memory requirements. For the forward mode the core memory 
increase of the AD-enhanced code is approximately NDV times the 
core memory requirement of the original (undifferentiated) code. 
For the reverse mode the corresponding increase is NOF time that 
of the original code. 

B. Second-Order Sensitivity Derivatives, 

SO Method 3, Airfoil Example 

Results are presented subsequently from the successful verifica- 
tion of the proposed efficient noniterative SO method 3 [Eq. (22)] 
for computing SO SDs. The example problem is steady transonic in- 
viscid flow over a NACA 0012 airfoil with fireestream Mach number 
Afoo 0.80 and angle of attack: a 1 deg. The two-dimensional Euler 
equations are solved on a Sun workstation in double precision using 
a conventional finite volume upwind flux- vector-splitting scheme. A 
C-mesh computational grid is used with dimensions 129 x 33 grid 
points. High-quality lift-corrected boundary conditions are used at 
the far field boundary, which is placed approximately five chord 
lengths from the surface of the airfoil. 

In the present example derivatives of three aerodynamic out- 
put functions are considered: C*., Cp, and Cm (Le., coefficients of 
lift, wave drag, and pitching moment, respectively). The computed 
steady-state values of these aerodynamic force coefficients are given 
in Tkble 2. Note that the number of digits shown in Table 2 (and also 
Tables 3 and 4) is to illustrate consistency rather than accuracy. In 
addition, derivatives with respect to three aerodynamic input vari- 
ables are considered. They arc g (a geometric-shape variable), a , and 
Afoo. The geometric-shape variable g is a single arbitrarily selected 
y coordinate of the computational grid on the surface of the airfoil: 
simple, but one that is sufficient fen* verifying geometric-shape SDs. 

Calculation of SO SDs by SO method 3 requires that all FO 
SDs are calculated first using both the forward-mode [Eqs. (3) and 
(4)] and the reverse-mode [Eqs. (7) and (8)] approaches. The cal- 
culated FO SDs from a hand-differentiated incremental-iterative 
(HDII) implementation of these two approaches are presented in 
Table 3, where the results are seen to agree, as expected. The FO 
SDs presented in Table 3 have been thoroughly verified for consis- 
tency through a meticulous implementation of the method of central 
finite differences, where agreement to six significant digits or greater 
is noted in all comparisons. 


Table 2 Aerodynamic force 
coefficients, airfoil example 


Coefficient 

Value 

Cl 

+0.2830659£+00 

C d 

+0.2070493E— 01 

Cm 

— 0.2876639E— 01 
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Ifcble3 First-order s e nsitivit y derivatives, 
HD1I method, airfoil example 


Derivative 

bj 

Forward mode 

Reverse mode 

dC L 

8 

+0.14G5406£+00 

+0.1405406£+00 

dbi 

a 

— 0.1087323£— 01 

—0. 1087323F— 01 

° °i 

Moo 

— 0.4672729£— 01 

— 0.4672729£ —01 

dCp 

g 

+0.1761807£+02 

+0.1761807£+02 

dbi 

a 

+0.1 158625F+01 

+0.1158625F+01 

j 

Mqo 

— 033825 80£+01 

— 03382580£ +01 

dC&t 

8 

+0.3171492£+01 

+03171492F+01 

db7 

a 

+03955598F+O1 

+03955598£ +01 

Moo 

— 0.1603002F+01 

-0.1603002F+01 


Mle4 Second-order snsttvity derivatives, 
SO method 3, airfoti example 


Derivative bj/bt g a M*, 

d 2 ^ g +0.248807F+Q3 +0.250402F+03 +0.2G5825F+03 

j! jl a +0.25G4G2E+03 +G.i84277£+05 +0.160858£+05 

J A#*, +Q-205825£+03 +0.160858£+05 +0.133087£+05 

<FCd g +0.71777£+02 +0.134379£+02 +0.101304£+02 

J, » a +0.134379F+02 +0.959310£+03 +0.804021 £+03 

/ 4 Moo +0.101304£+02 +0.804021 £+03 +0.662088£+O3 

&Ch 8 — 0.663590£ +02 -0.602441 £+02 -0.490399£+02 

a -0.6024415+02 -0.449198£+04 -0386943+04 
J k Moo — 0.490399£ +02 -0386943£+04 -03205 12£+04 


Table 5 Relative CPU timings: complete SO 
method 3* airfoil example 


Computational procedure 

% of total 

Nonlinear flow, Eqs. (1) and (2) 

5.4 

Forward-mode FO SDs, Eqs. (3) and (4) 

253 

Reverse-mode FO SDs, Eqs. (7) and (9) 

69.0 

SO SDs, Eq. (22) 

0.1 

Total 

100.0 


Table 5 shows that using an ADIFOR-assisted second differentia- 
tion SO SDs can be obtained extremely fast, if one already has both 
the forward-mode and reverse-mode FO SDs. 


IV. Conclusions 

An efficient incremental-iterative approach for differentiating ad- 
vanced CFD flow codes has been successfully demonstrated on a 
two-dimensional inviscid model problem. The method employs the 
reverse-mode capability of the automatic-differentiation software 
tool ADIFOR 3.0 and has been shown to yield consistent first-order 
aerodynamic sensitivity derivatives. A substantial reduction in CPU 
time and computer memory has been demonstrated by comparison 
with results from a straightforward, black-box reverse-mode appli- 
cation of ADIFOR 3.0 to die same flow code. 

A computationally efficient ADIFQR-assisted procedure for con- 
sistent second-order aerodynamic sensitivity derivatives has been 
successfully verified on an inviscid transonic lifting airfoil exam- 
ple problem. Accurate second derivatives (i.e., the complete Hes- 
sian matrices) of lift, wave drag, and pitching-moment coefficients 
with respect to geometric shape, angle of attack, and ffeestream 
Mach number have been calculated. With the present procedure 
second-order derivatives are now computationally feasible, at least 
in two dimensions. The computation of secrakU>rder derivatives 
in three dimensions appeals to be within reach, but remains to be 
investigated. 

This second-order method requires that first-order derivatives be 
calculated using both the forward (direct) and reverse (adjoint) pro- 
cedures; then, second-order derivatives can be obtained in a noniter- 
ative calculation that is computationally very efficient. An ADIFOR 
differentiation is used to generate a number of required second-order 
terms in this noniterative calculation. If one already has either for- 
ward (NDV solutions) or reverse (NOF solutions) FO SDs, then 
upon obtaining the other FO SDs (NOF or NDV additional solu- 
tions, respectively) one calculates all of die SO SDs ( NOF x NDV 2 
derivatives) very efficiently. 


The SO method 3 is implemented by application (in the for- 
ward mode) of ADIFOR to appropriate pieces of the FORTRAN 
code used earlier for hand-differentiated forward-mode calculation 
of the FO SDs. The calculated SO SDs from this implementation 
of SO method 3 are presented in Thble 4. The SO SDs of Table 
4 have been thoroughly verified for consistency through a meticu- 
lous application of central finite differences applied to FO SDs ob- 
tained by the hand-differentiated methods already described. Agree- 
ment to five significant digits or better is noted in this verification 
study for all SO SDs reported in Table 4. This verification study 
was not conducted using finite differences applied to the origi- 
nal nonlinear flow code; that approach has been documented to 
be vulnerable to severe numerical inaccuracy when SO SDs are 
calculated. 3 The symmetry of the calculated SO SDs shown in Th- 
ble 4 is expected and results from the computations perf orm ed, 
that is, no derivative symmetry was explicitly imposed on the 
problem. 

For the present airfoil example problem Table 5 illustrates (in 
terms of percentages of the total) the breakdown of relative CPU 
timings for the computational steps of SO method 3 procedure for 
calculating the SO SDs. Not included in Thble 5 is the CPU time 
for the grid generation and the grid-sensitivity derivatives, negligi- 
ble for this particular two-dimensional example. Thble 5 illustrates 
clearly the computational efficiency of the SO method 3 for SO SDs. 
Recall that results of the present example are for three aerodynamic 
output functions and three input (design) variables, where the com- 
putational work of the forward-mode and reverse-mode procedures 
for FO SDs should be approximately equal (in theory, for hand- 
differentiated code, as used here). In this example, however. Table 5 
reveals that the reverse mode was much more costly than the for- 
ward mode; apparently the three linear systems for the reverse mode 
are suffer than the three for the forward mode. This characteristic 
of the adjoint equations has been observed by others. 7 As expected. 


Appendix: Expansion of Terms 

In this Appendix the terms DF T y /Dfr*, DFJ./Db*, and DG, m /D&t 
are expanded using the index notation alrcaay established. The ex- 
pansion of VF'j/Dbt is 


Vb k 


= — Q'„ 

8Q, " 


nk 


A' 

+ ax, 





a 2 F, 

3G.ae« 


Q'm,Q'* + 


8 2 F i 

bQ,8X p 


K,Q'« 


+ 


d 2 Fi 

dQndbj 


Q'm 


, th O' X' , X' X' +3!£Lx 

+ ax,a Q„ mi 9 k+ dx 9 dx p ri qk + dx q dbj qk 


+ 


a 2 Fi 

*b k 8Q m 


Q-J + 


#Fj , 8 2 Fi 

8b k 8X r pi db t dbj 


(Al) 


In Eq. (Al), the indices i, j, and k are free and repeated indices it, 
m, p, and q are summed. The terms of DRy/Dbt are obtained from 
Eq. (Al) by replacing everywhere F'j with R[j and F, with Ri (and 
thus / replaces i as a free index in the resulting expressions). 
Finally, the expansion of the terms for DG, n /Dfr t is 


DG im aG iW . BG im dG in 

Db t dQ, Unk ax„ 8b k 


_J!^ eu+ _i2L_ 

aQ,8Qm * bQ,8Q m 


Q'« + hi 


9 2 *< y) 
8X q 8Q m «* 


, *F t 8 2 R, d 2 F t 

BX q 8Q m * k "db k 8Q m 8b k 8Q m 


(A2) 
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In Eq. (A2) the indices i, m, and k are free, and repeated indices /, 
ft, and q are summed. 
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